module mesh_manager

    use phy_state
    use interp
    use mesh

    implicit none

    private
    
    public &
        mesh_manager_location_query, &
        mesh_manager_boundary_check, &
        mesh_manager_variable_query, &
        mesh_manager_divergence_query

contains

    subroutine mesh_manager_location_query(x, ci)
        real, intent(in) :: x(num_dim)
        integer, intent(out) :: ci

        integer i, o(num_dim)

        if(is_domain_regular) then
            if(is_mesh_uniform) then
                do i = 1, num_dim
                    o(i) = int(x(i)/dx(i))+1
                    if(o(i) > num_cell_(i)) o(i) = o(i)-1
                end do
                ci = cell_index_stretch(o)
            end if
        else
            ! ---> generic mesh to be constructed
        end if

        return

    contains

        integer function cell_index_stretch(o) result(res)
            integer, intent(in) :: o(num_dim)
            
            integer i
            
            res = o(1)
            do i = 2, num_dim
                res = res+(o(i)-1)*num_slice_cell(i)
            end do
            
            return
        end function cell_index_stretch

    end subroutine mesh_manager_location_query

    subroutine mesh_manager_boundary_check(x)
        real, intent(inout) :: x(num_dim)

        integer i

        if(is_domain_regular) then
            do i = 1, num_dim
                if(x(i) < boundary_x(1,i)) then
                    x(i) = x(i)+span(i)
                else if(x(i) > boundary_x(2,i)) then
                    x(i) = x(i)-span(i)
                end if
            end do
        else
            ! ---> generic mesh to be constructed
        end if

        return
    end subroutine mesh_manager_boundary_check

    subroutine mesh_manager_variable_query(var_name, it, x, v, ci)
        character(*), intent(in) :: var_name
        integer, intent(in) :: it
        real, intent(in) :: x(num_dim)
        real, intent(out) :: v
        integer, intent(in), optional :: ci

        real, allocatable :: vertex_v(:)
        integer i, j, k

        if(present(ci)) then
            i = ci
        else
            call mesh_manager_location_query(x, i)
        end if

        allocate(vertex_v(cells(i)%num_vertex))
        do j = 1, cells(i)%num_vertex
            k = cells(i)%vertex_index(j)
            vertex_v(j) = grids(k)%pstate%variable_query(var_name, it)
        end do

        call interp_linear(num_dim, cells(i)%num_vertex, &
            cells(i)%vertex_x, vertex_v, x, v)

        return
    end subroutine mesh_manager_variable_query

    subroutine mesh_manager_divergence_query(it, x, div, ci)
        integer, intent(in) :: it
        real, intent(in) :: x(num_dim)
        real, intent(out) :: div
        real, intent(in), optional :: ci

        real, allocatable :: vertex_v(:,:)
        real dudx, dvdy, dwdz
        integer i, j, k

        if(present(ci)) then
            i = ci
        else
            call mesh_manager_location_query(x, i)
        end if

        allocate(vertex_v(num_dim,cells(i)%num_vertex))
        do j = 1, cells(i)%num_vertex
            k = cells(i)%vertex_index(j)
            if(num_dim >= 1) then
                vertex_v(1,j) = grids(k)%pstate%variable_query("u", it)
            end if
            if(num_dim >= 2) then
                vertex_v(2,j) = grids(k)%pstate%variable_query("v", it)
            end if
            if(num_dim >= 3) then
                vertex_v(3,j) = grids(k)%pstate%variable_query("w", it)
            end if
        end do

        div = 0.0d0
        if(num_dim >= 1) then
            call interp_linear_derv_x(num_dim, cells(i)%num_vertex, &
                cells(i)%vertex_x, vertex_v(1,:), x, dudx)
            div = div+dudx
        end if
        if(num_dim >= 2) then
            call interp_linear_derv_y(num_dim, cells(i)%num_vertex, &
                cells(i)%vertex_x, vertex_v(2,:), x, dvdy)
            div = div+dvdy
        end if
        if(num_dim >= 3) then
            call interp_linear_derv_z(num_dim, cells(i)%num_vertex, &
                cells(i)%vertex_x, vertex_v(3,:), x, dwdz)
            div = div+dwdz
        end if

        return
    end subroutine mesh_manager_divergence_query

end module mesh_manager
